QC using RSEM data

counts <- read_csv("data/genes_counts.csv")
tpm <- read_csv("data/genes_tpm.csv")
treatments <- read_csv("data/treatments.csv", col_types = "ccc")
ercc <- left_join(read_tsv("data/ercc.tsv"), tpm, by = "gene_id") 
descriptions <- read_tsv("data/descriptions.tsv", col_names = c("gene_id", "description"))
mix1 <-  ercc[,c("gene_id", "mix1", subset(treatments, spikein == "mix1")$sample)] %>% gather( library, observed, -c(mix1,gene_id)) %>% dplyr::rename(expected = mix1 )
mix2 <- ercc[,c("gene_id", "mix2", as.character(subset(treatments, spikein == "mix2")$sample))] %>% gather( library, observed, -mix2, -gene_id) %>% dplyr::rename(expected = mix2 )
ercc <- rbind(mix1,mix2)
ggplot(ercc, aes(observed+1, expected+1, color=as.factor(library))) + geom_point() + scale_y_continuous(trans='log2')+scale_x_continuous(trans='log2')+ geom_abline()+guides(color=F)

cutoff <- seq(0,5,by=0.1)
r <- c()
for (i in cutoff) {
  keepGenes <- ercc %>% group_by(gene_id, library) %>% dplyr::summarise(all_mean = mean(observed)) %>% dplyr::filter(all_mean > i) %>% .$gene_id %>% unique
  r <- c(r, summary(lm(log2(observed+1)~log2(expected+1), data=subset(ercc, gene_id %in% keepGenes)))$r.squared)
}
plot(cutoff,r)

The signal to noise is basically fine even at the lowest levels, so we do not do any abundance-based filtering.

Table S1: Correlation between pheromone components and principal component axes

loadings <- pca.all$rotation
sweep(x = abs(loadings), MARGIN = 2, STATS = colSums(abs(loadings)), FUN = "/")[,1:2] %>% kable(output="html", digits=3)

Figure S2: Correlation matrix between individual pheromone components

M<-cor(qp[, 1:6])

colnames(M) <- rownames(M) <- c( "10-HDAA", "9-HDA", "HOB", "10-HDA", "HVA", "9-ODA")

cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

p.mat <- cor.mtest(qp[, 1:6]) # p-value correlation matrix

col <- colorRampPalette(c("#4477AA", "#77AADD", "#FFFFFF", "#EE9988","#BB4444"))

corrplot::corrplot(M, method="color", col=col(200),  
         type="lower", order="hclust", outline=T, mar = c(1, 1, 1, 7),
         addCoef.col = "black", tl.col="black", tl.srt=80, p.mat = p.mat, sig.level = 0.01, insig = "n",diag=FALSE)

Choose the number of principal components for principal component regression

first <- clmm2(ovary ~ Treatment*(PC1), data = qp2)
second <- clmm2(ovary ~ Treatment*(PC1+PC2), data = qp2)
third <- clmm2(ovary ~ Treatment*(PC1+PC2+PC3), data = qp2)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

anova(first,second)

## Likelihood ratio tests of cumulative link models
## 
## Response: ovary
##                           Model Resid. df -2logLik   Test    Df LR stat.
## 1       Treatment * (PC1) |  |         28 79.15867                      
## 2 Treatment * (PC1 + PC2) |  |         26 70.23393 1 vs 2     2 8.924741
##      Pr(Chi)
## 1           
## 2 0.01153499

anova(second,third)

## Likelihood ratio tests of cumulative link models
## 
## Response: ovary
##                                 Model Resid. df -2logLik   Test    Df
## 1       Treatment * (PC1 + PC2) |  |         26 70.23393             
## 2 Treatment * (PC1 + PC2 + PC3) |  |         24 68.41331 1 vs 2     2
##   LR stat.   Pr(Chi)
## 1                   
## 2 1.820621 0.4023992

Adding the first two principal components improves overall model fit, adding the third does not

ggplot(qp2, aes(as.factor(`Ovary stage`), PC2, fill=Treatment))+geom_line(stat="smooth", linetype = 2, method = "rlm", se=F, aes(color= Treatment, group=Treatment, alpha=.1))+geom_boxplot(outlier.color="white", varwidth = TRUE, position=position_dodge2(preserve = "single"))+theme_bw()+ylab("Queen mandibular pheromone (PC2)")+xlab("Ovary activation stage")+scale_fill_manual(values=c("steelblue3", "orange"))+scale_color_manual(values=c("steelblue3", "orange"))+guides(alpha = F)+theme(legend.position="top", legend.text=element_text(size=12), axis.text=element_text(size=14))+coord_trans(limy=c(-1.5,.7))+theme(axis.title=element_text(size=14, face="bold"))

#ggsave("figures/qp ovary.pdf", height=5, width=5) # raw Figure 1

Table S2: Cumulative link mixed model analysis of pheromonal composition and ovary development

We see that treatment has a strongly significant effect on ovary activation, and there is an interaction between PC2 and treatment and their relationship to ovary development

clmm.mod <- clmm2(ovary ~ Treatment*(PC1+PC2), data = qp2)
summary(clmm.mod)

## Call:
## clm2(location = ovary ~ Treatment * (PC1 + PC2), data = qp2)
## 
## Location coefficients:
##                  Estimate Std. Error z value Pr(>|z|) 
## TreatmentGFP      2.8107   0.9148     3.0724 0.0021232
## PC1              -0.3679   0.4404    -0.8353 0.4035230
## PC2              -2.6820   1.2409    -2.1614 0.0306667
## TreatmentGFP:PC1  0.1099   0.4972     0.2210 0.8251167
## TreatmentGFP:PC2  3.9023   1.6861     2.3144 0.0206461
## 
## No scale coefficients
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 0|1  0.6350   0.5807     1.0935
## 1|2  1.7584   0.6685     2.6302
## 2|3  2.5175   0.7378     3.4120
## 
## log-likelihood: -35.11696 
## AIC: 86.23393 
## Condition number of Hessian: 240.3247

Table S3: There are no differences in mortality between the treatments

bm <- read_tsv("data/DSX_RNAi_Mortality.txt", col_types = cols())
summary(glm(Mortality ~ 1 + Treatment,  data = bm, binomial(link= "logit")))

## 
## Call:
## glm(formula = Mortality ~ 1 + Treatment, family = binomial(link = "logit"), 
##     data = bm)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.328  -1.310   1.034   1.051   1.051  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.30538    0.35222   0.867    0.386
## TreatmentGFP  0.04293    0.51596   0.083    0.934
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 84.330  on 61  degrees of freedom
## Residual deviance: 84.323  on 60  degrees of freedom
## AIC: 88.323
## 
## Number of Fisher Scoring iterations: 4

Figure 2

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:Biobase':
## 
##     combine

## The following object is masked from 'package:BiocGenerics':
## 
##     combine

## The following object is masked from 'package:dplyr':
## 
##     combine

ovaryPlot <- ggplot(qp2, aes(x=as.factor(`Ovary stage`), fill=Treatment))+geom_histogram(stat="count",position=position_dodge2(preserve = "single"))+scale_fill_manual(values=c("steelblue3", "orange"))+theme_bw()+coord_flip()+xlab("Ovary activation stage")

## Warning: Ignoring unknown parameters: binwidth, bins, pad

qpPlot <- ggplot(qp2, aes(as.factor(`Ovary stage`), PC2, fill=Treatment))+geom_line(stat="smooth", linetype = 2, method = "rlm", se=F, aes(color= Treatment, group=Treatment, alpha=.1))+geom_boxplot(outlier.color="white",  position=position_dodge2(preserve = "single"))+theme_bw()+ylab("Queen mandibular pheromone (PCA2)")+xlab("")+scale_fill_manual(values=c("steelblue3", "orange"))+scale_color_manual(values=c("steelblue3", "orange"))+guides(alpha = F)+coord_flip(ylim=c(-1.5,.7))

pheroMeans <- qp %>% dplyr::select(`10-HDAA`, `9-HDA`,HOB,`10-HDA`,HVA,`9-ODA`)  %>% colMeans()
pheroMeans <- data.frame(mean = pheroMeans, compound = names(pheroMeans))

loadingsPlot <-   sweep(x = abs(loadings), MARGIN = 2, STATS = colSums(abs(loadings)), FUN = "/")[,1:2] %>%  as.data.frame %>% rownames_to_column(var="compound") %>% gather(key="PC", value="loadings", -compound) %>% mutate(type = ifelse(compound %in% c("HOB", "9-HDA", "HVA", "9-ODA"), "queen", "worker")) %>% left_join(pheroMeans, by="compound")  %>% filter(PC == "PC2") %>% ggplot(aes(PC,compound))+geom_point(aes(fill=loadings, color=type,size=mean),pch=21,stroke = 1)+scale_fill_gradient(low="white",high="magenta")+scale_color_manual(values=c("red","blue"))+theme_bw()+theme(axis.title.y=element_blank()) + scale_size_continuous(range = c(2,6))+guides(size=F)+ coord_flip()+ theme(legend.direction="horizontal")

## Warning: Column `compound` joining character vector and factor, coercing
## into character vector

pdf("figures/Fig2.pdf", height=5,width=7, useDingbats=F)
grid.arrange(arrangeGrob(ggarrange(ovaryPlot, qpPlot, labels = c("A", "C"), ncol = 2, nrow = 1, common.legend = T )), arrangeGrob(loadingsPlot), ncol=1, heights=c(3, 1))
dev.off()

## quartz_off_screen 
##                 2

Table S4: Distribution of genes by module. Dsx containing module is highlighted in red

Exploring the Dsx-sensitive module

Alldegrees <- intramodularConnectivity.fromExpr(expDat, net$colors, power = SoftPower)

##  softConnectivity: FYI: connecitivty of genes with less than 7 valid samples will be returned as NA.
##  ..calculating connectivities..

rownames(Alldegrees) <- colnames(expDat)
geneModuleMembership <- as.data.frame(cor(expDat, MEs, use = "p"));
Alldegrees["Dsx",]

##       kTotal  kWithin     kOut     kDiff
## Dsx 76.54159 32.44871 44.09288 -11.64418

DsxModuleGenes <- colnames(expDat)[net$colors == DsxModule]
connectivityDsxModule <- data.frame(kWithin = Alldegrees[DsxModuleGenes,"kWithin"], gene_id = rownames(Alldegrees[DsxModuleGenes,]), moduleMembership = geneModuleMembership[DsxModuleGenes, paste0("ME",DsxModule)]) %>% left_join(de.tgw$table, by = "gene_id")

## Warning: Column `gene_id` joining factor and character vector, coercing
## into character vector

dsxModuleGenes <- descriptions %>% filter(gene_id %in% colnames(expDat)[net$colors == DsxModule] ) %>% left_join(de.tgw$table, by = "gene_id") %>% left_join(connectivityDsxModule %>% dplyr::select(gene_id, kWithin), by="gene_id") %>% dplyr::select(gene = gene_id, description = description.x, logFC, kWithin, logCPM) %>% arrange(desc(logFC))  
highligh <- which(dsxModuleGenes$gene %in% c("Vg", "Dsx", "Dnmt3"))